(* ::Package:: *)

(*******************************************************************
This file was generated automatically by the Mathematica front end.
It contains Initialization cells from a Notebook file, which
typically will have the same name as this file except ending in
".nb" instead of ".m".

This file is intended to be loaded into the Mathematica kernel using
the package loading commands Get or Needs.  Doing so is equivalent
to using the Evaluate Initialization Cells menu command in the front
end.

DO NOT EDIT THIS FILE.  This entire file is regenerated
automatically each time the parent Notebook file is saved in the
Mathematica front end.  Any changes you make to this file will be
overwritten.
***********************************************************************)















BeginPackage["MathSVM`"];
(*"Graphics`Graphics`",{"LinearAlgebra`MatrixManipulation`"}*)


SVMClassify::usage="Classify[K,X,\[Alpha],y,x] predicts a response y for a new case x using the trained svm model {kernel K, data X,multipliers \[Alpha],labels y}.";
KernelMatrix::usage="KernelMatrix[K,X] computes the (full) matrix with elements K[x\[LeftDoubleBracket]i\[RightDoubleBracket],x\[LeftDoubleBracket]j\[RightDoubleBracket]].";
IdentityKernel::usage="IdentityKernel[x,y] is the function x.y.";
RBFKernel::usage="RBFKernel[x,y] is the function Exp[-\[Gamma](x-y).(x-y)].";
PolynomialKernel::usage="PolynomialKernel[x,y,d], where d is an integer, is the function (x.y+1\!\(\*SuperscriptBox[\")\", \"d\"]\).";
SeparableSVM::usage="SeparableSVM[X,y,\[Tau]] trains a separable SVM on data X, labels y, and solution tolerance \[Tau] (see QPSolve). Returns the multiplier vector \[Alpha]. Option KernelFunction determines the kernel to use; default is IdentityKernel (i.e., no kernel).";
NonseparableSVM::usage="NonSeparableSVM[X,y,C,\[Tau]] trains an non-separable SVM on data X, labels y, penalty term C and solution tolerance \[Tau] (see QPSolve). Returns the multiplier vector \[Alpha]. Option KernelFunction determines the kernel to use; default is IdentityKernel (i.e., no kernel).";
RegressionSVM::usage="RegressionSVM[X,y,C,\[Epsilon],\[Tau]] trains a regression SVM with \[Epsilon]-insensitive loss function on regression variables X, response variable y, error tolerace C and solution tolerance \[Tau] (see QPSolve). Returns the multiplier vector \[Alpha]. Option KernelFunction determines the kernel to use; default is IdentityKernel (i.e., no kernel). Note that this creates a dual problem of size 2*(#samples).";
RegressionBias::usage="RegressionBias[\[Alpha],X] returns the bias term for regression problems (this differs form the classification formulations).";
RegressionFunction::usage="RegressionFunction[\[Alpha],X,x] returns the regression function f(x) defined by \[Alpha] and X. Option KernelFunction determines the kernel to use; default is IdentityKernel (i.e., no kernel).";
RegressionSupportVectors::usage="RegressionSupportVectors[\[Alpha]] returns the support vectors for a regression model defined by \[Alpha].";
SupportVectors::usage="SupportVectors[\[Alpha],y] returns {I+,I-} where I+/- are the index vector for support vectors in class +/-.";
WeightVector::usage="WeightVector[\[Alpha],X,y] returns the weight (normal) vector for the separating plane. This is only well-defined when \[Alpha] is the result of an SVM trained with IdentityKernel.";
Bias::usage="Bias[\[Alpha],X,y] gives the bias term for the trained svm model {kernel K, data X,multipliers \[Alpha],labels y}. Option KernelFunction determines the kernel to use; default is IdentityKernel (i.e., no kernel).";
KernelDistance::usage="KernelDistance[x,y,K] gives the euclidean (\!\(\*SubscriptBox[\"L\", \"2\"]\)) distance between x and y in the feature space induces by the kernel K.";
SVMPlot::usage="SVMPlot[\[Alpha],X,y] produces a typical SVM-plot of the SVM solution \[Alpha], together with data points X (in \!\(\*SuperscriptBox[\"R\", \"2\"]\)!) and labels y. Option KernelFunction determines the kernel to use for drawing decision boundaries.";
SVMDataPlot::usage="SVMDataPlot[X,y] just plots the data points in X with classes {+1,-1} highlighted.";
RegressionSVMPlot::usage="RegressionSVMPlot[\[Alpha],X,y,\[Epsilon]] draws data points {X,y}, the regression function determined by {\[Alpha],X} and a 'corridor' of with \[Epsilon].";
QPSolve::usage="QPSolve[Q,p,a,b,c,y,\[Tau]] solves the quadratic programming problem min \[Alpha].Q.\[Alpha]+p.\[Alpha], subject to a\[LessEqual]\[Alpha]\[LessEqual]b and y.\[Alpha]=c. QPSolve uses the GSMO algorithm described by Keerthi et al. \[Tau] is a solution tolerance parameter (0.01 or so is usually good enough for SVMs). Q must be a positive semidefinite matrix to guarantee convergence.";



FormulationOnly::usage="This option causes SVM training functions to merely formulate the QP, without solving it.";\

KernelFunction::usage="KernelFunction->K[x,y] is used in many SVM-related functions to define the kernel function.";\

DebugSolver;

Begin["MathSVM`Private`"];





KernelMatrix[K_,X_]:=Outer[K,X,X,1]



IdentityKernel[x_,y_]=x.y







RBFKernel[x_,y_,\[Gamma]_]=Exp[-\[Gamma](x-y).(x-y)]





PolynomialKernel[x_,y_,d_Integer]=(x.y+1)^d







SVMClassify[K_,X_,\[Alpha]_,y_,b_,x_]:=
  Sign[Total[y*\[Alpha]*Map[K[#,x]&,X]]+b]

SVMClassify[K_,X_,\[Alpha]_,y_,x_]:=
  Sign[Total[y*\[Alpha]*Map[K[#,x]&,X]]+Bias[\[Alpha],X,y,K]]





Options[SeparableSVM]={FormulationOnly->False,
    KernelFunction->IdentityKernel}

SeparableSVM[X_,y_,\[Tau]_,opts___]:=
  Module[{l,qp,K},
    l=Length[X];
    K=KernelFunction/.{opts}/.Options[SeparableSVM];
    qp=
      {KernelMatrix[K,X]*Outer[Times,y,y],
        Table[-1,{l}],
        Table[0,{l}],
        Table[\[Infinity],{l}],
        0,
        y,
        \[Tau]};
    If[FormulationOnly/.{opts}/.Options[SeparableSVM],
      qp,
      Apply[QPSolve,qp]]]







Options[NonseparableSVM]={FormulationOnly->False,
    KernelFunction->IdentityKernel}

NonseparableSVM[X_,y_,C_,\[Tau]_,opts___]:=
  Module[{l,qp,K},
    l=Length[X];
    K=KernelFunction/.{opts}/.Options[NonseparableSVM];
    qp=
      {KernelMatrix[K,X]*Outer[Times,y,y],
        Table[-1,{l}],
        Table[0,{l}],
        Table[C,{l}],
        0,
        y,
        \[Tau]};
    If[FormulationOnly/.{opts}/.Options[NonseparableSVM],
      qp,
      Apply[QPSolve,qp]]]











Options[RegressionSVM]={FormulationOnly->False,
    KernelFunction->IdentityKernel}

RegressionSVM[X_,z_,C_,\[Epsilon]_,\[Tau]_,opts___]:=
  Module[{l,Q,\[Alpha]},
    l=Length[X];
    K=KernelFunction/.{opts}/.Options[RegressionSVM];
    qp=
      {KernelMatrix[K,Join[X,-X]],
        Table[\[Epsilon],{2l}]+Join[z,-z],
        Join[Table[0,{2l}]],
        Join[Table[C,{2l}]],
        0,
        Join[Table[1,{l}],Table[-1,{l}]],
        \[Tau]};
    If[FormulationOnly/.{opts}/.Options[RegressionSVM],
      qp,
      Apply[QPSolve,qp]]]

RegressionSVM[X_?VectorQ,z_,C_,\[Epsilon]_,\[Tau]_,opts___]:=
  RegressionSVM[Partition[X,1],z,C,\[Epsilon],\[Tau],opts]



RegressionBias[\[Alpha]_,X_?MatrixQ,z_,\[Epsilon]_,opts___]:=
  Module[{l=Length[X],K},
    K=KernelFunction/.{opts}/.Options[RegressionSVM];
    sv=SupportVectors[\[Alpha]][[1]];
    \[Epsilon]+z[[sv]]-
      Sum[(\[Alpha][[
                i+l]]-\[Alpha][[
                i]])K[
            X[[i]],
            X[[sv]]],{i,l}]]

RegressionBias[\[Alpha]_,X_?VectorQ,z_,\[Epsilon]_,opts___]:=
  RegressionBias[\[Alpha],Partition[X,1],z,\[Epsilon],opts]

RegressionFunction[\[Alpha]_,X_?MatrixQ,z_,\[Epsilon]_,x_?VectorQ,opts___]:=
  Module[{l=Length[X],K,b},
    b=RegressionBias[\[Alpha],X,z,\[Epsilon],opts];
    K=KernelFunction/.{opts}/.Options[RegressionSVM];
    Sum[(\[Alpha][[
                i+l]]-\[Alpha][[
                i]])K[
            X[[i]],x],{i,l}]+b]

RegressionFunction[\[Alpha]_,X_?VectorQ,z_,\[Epsilon]_,x_,opts___]:=
  RegressionFunction[\[Alpha],Partition[X,1],z,\[Epsilon],{x},opts]



RegressionSupportVectors[\[Alpha]_List]:=
  Module[{l=Length[\[Alpha]]/2},
    SupportVectors[\[Alpha],Join[Table[1,{l}],Table[-1,{l}]]]-{0,l}]







SupportVectors[\[Alpha]_List,
    y_List]:={Position[\[Alpha]*y,_?Positive]//Flatten,
    Position[\[Alpha]*y,_?Negative]//Flatten}



SupportVectors[\[Alpha]_List]:=
  Flatten[Position[\[Alpha],a_/;a!=0]]



WeightVector[\[Alpha]_List,X_?MatrixQ,y_List]:=
  Sum[y[[i]]*\[Alpha][[
        i]]*X[[i]],{i,
      Length[X]}]



Options[Bias]={KernelFunction->IdentityKernel}

Bias[\[Alpha]_,X_,y_,opts___]:=
  Module[{sv,K},
    K=KernelFunction/.{opts}/.Options[Bias];
    sv=SupportVectors[\[Alpha],y];
    -(Total[\[Alpha]*y*
                Map[K[X[[
                        sv[[1,
                          1]]]],#]&,X]]+
            Total[\[Alpha]*y*
                Map[K[X[[
                        sv[[2,
                          1]]]],#]&,
                  X]])/2]





KernelDistance[x_,y_,K_]:=Sqrt[K[x,x]+K[y,y]-2K[x,y]]







Options[SVMPlot]={KernelFunction->IdentityKernel}

SVMPlot[\[Alpha]_,X_?MatrixQ,y_,opts___]:=
  Module[{x1range,x2range,sv,df,x1,x2,K},
    K=KernelFunction/.{opts}/.Options[SVMPlot];
    x1range={x1,Min[X[[All,1]]],
        Max[X[[All,1]]]};
    x2range={x2,Min[X[[All,2]]],
        Max[X[[All,2]]]};
    sv=SupportVectors[\[Alpha],y];
    df=Total[\[Alpha]*y*Map[K[{x1,x2},#]&,X]]+Bias[\[Alpha],X,y,opts];
    DisplayTogether[
      ListPlot[X[[Flatten[sv]]],
        PlotStyle->{Hue[0.6],PointSize[0.015]}],
      ListPlot[Extract[X,Position[y,1]]],
      ListPlot[Extract[X,Position[y,-1]],PlotStyle->{GrayLevel[0.6]}],
      ContourPlot[df==0,x1range,x2range],
      ContourPlot[df==1,x1range,x2range,
        PlotStyle->Dashing[{.01,.01}]],
      ContourPlot[df==-1,x1range,x2range,
        PlotStyle->Dashing[{.01,.01}]],
      RemainingOptions[{KernelFunction},opts]]]

SVMDataPlot[X_?MatrixQ,y_,opts___]:=
  DisplayTogether[
    ListPlot[Extract[X,Position[y,1]]],
    ListPlot[Extract[X,Position[y,-1]],PlotStyle->{GrayLevel[0.6]}],opts]



Options[RegressionSVMPlot]={KernelFunction->IdentityKernel}

RegressionSVMPlot[\[Alpha]_,X_?VectorQ,z_,\[Epsilon]_,opts___]:=
  Module[{K,rf,sv},
    K=KernelFunction/.{opts}/.Options[RegressionSVMPlot];
    rf=RegressionFunction[\[Alpha],X,z,\[Epsilon],x,opts];
    sv=RegressionSupportVectors[\[Alpha]];
    DisplayTogether[
      Plot[rf,{x,Min[X],Max[X]}],
      Plot[rf+\[Epsilon],{x,Min[X],Max[X]},
        PlotStyle->Dashing[{.01,.01}]],
      Plot[rf-\[Epsilon],{x,Min[X],Max[X]},
        PlotStyle->Dashing[{.01,.01}]],
      ListPlot[
        Thread[{X[[Flatten[sv]]],
            z[[Flatten[sv]]]}],
        PlotStyle->{Hue[0.6],PointSize[0.015]}],
      ListPlot[Thread[{X,z}]],
      PlotRange->All]]











Options[QPSolve]={DebugSolver->False}

QPSolve[Q_,p_,a_,b_,c_,y_,\[Tau]_,opts___]:=
  Module[{\[Alpha],l,Isets,Iup,Ilow,F,B,oldB,M,k,debug},
    debug=DebugSolver/.{opts}/.Options[QPSolve];
    l=Length[Q];
    \[Alpha]=FeasiblePoint[a,b,y,c];
    If[\[Alpha]==Null,Return[Null]];k=0;
    While[True,
      Isets=GetIndexSets[\[Alpha],a,b,y];
      Iup=UpperBoundarySet[Isets];
      Ilow=LowerBoundarySet[Isets];
      F=(Q.\[Alpha]+p)/y;
      If[QPOptimalQ[Iup,Ilow,F,\[Tau]],Break[]];
      B=GetViolatingPair[Isets[[1]],Iup,
          Ilow,F,\[Tau]];
      If[B==Null,Return["Error: no violating pair"]];
      If[B==oldB,
        Print["Error: stuck on violating pair ",B];
        Print[
          "Subproblem = ",{Q[[B,B]],
            p[[B]]+
              Q[[B,
                  M]].\[Alpha][[
                  M]],
            a[[B]],
            b[[B]],
            y[[
              B]],\[Alpha][[
              B]],
            y[[
                B]].\[Alpha][[
                B]]}];
        Print["F = ",F[[B]]];
        Print["Isets = ", Isets];
        Return[\[Alpha]];];
      M=Complement[Range[1,l],B];
      If[debug,
        Print["--------------NEW PAIR:------------"];
        Print["viol. pair = ",B];
        Print[
          "\[Alpha] = ",\[Alpha][[B]]];
        Print["grad = ",F[[B]]];
        Print["current tau = ",Abs[
            F[[
                B[[1]]\
]]-
              F[[
                B[[2]]\
]]]]];
      \[Alpha][[B]]=
        QPSolve[Q[[B,B]],
          p[[B]]+
            Q[[B,
                M]].\[Alpha][[M\
]],a[[B]],
          b[[B]],
          y[[B]],\[Alpha]\
[[B]]];
      If[debug,
        Print["--------------NEW SOLUTION:------------"];
        Print[
          "\[Alpha] = ",\[Alpha][[B]]]];
      oldB=B];
    \[Alpha]]





FeasiblePoint[a_,b_,y_,c_]:=
  Module[{\[Alpha],l=Length[y]},
    \[Alpha]=Table[0,{l}];
    Do[
      \[Alpha][[
          i]]=(c-Drop[y,{i}].Drop[\[Alpha],{i}])/
          y[[i]];
      \[Alpha][[i]]=
        Min[Max[\[Alpha][[i]],
            a[[i]]],
          b[[i]]];
      If[y.\[Alpha]==c,Break],
      {i,l}];
    If[y.\[Alpha]==c,\[Alpha],Null]]





QPSolve[Q_/;Length[Q]==2,p_,a_,b_,y_,\[Alpha]_]:=
Module[{t,\[Alpha]new,c,tn,td},
c=\[Alpha].y;
tn=-((Q[[1]].\[Alpha]+p[[1]])/y[[1]]-(Q[[2]].\[Alpha]+p[[2]])/y[[2]]);
td=Q[[1,1]]/y[[1]]^2+Q[[2,2]]/y[[2]]^2-(Q[[1,2]]+Q[[2,1]])/(y[[1]]y[[2]]);
t=If[td!=0,tn/td,Sign[tn]*\[Infinity]];
\[Alpha]new=\[Alpha]+{t/y[[1]],-t/y[[2]]};
If[\[Alpha]new[[1]]<a[[1]],\[Alpha]new={a[[1]],(c-a[[1]]y[[1]])/y[[2]]}];
If[\[Alpha]new[[1]]>b[[1]],\[Alpha]new={b[[1]],(c-b[[1]]y[[1]])/y[[2]]}];
If[\[Alpha]new[[2]]<a[[2]],\[Alpha]new={(c-a[[2]]y[[2]])/y[[1]],a[[2]]}];
If[\[Alpha]new[[2]]>b[[2]],\[Alpha]new={(c-b[[2]]y[[2]])/y[[1]],b[[2]]}];
\[Alpha]new]





MyKronDelta[a_?VectorQ,b_?VectorQ]:=
  Map[KroneckerDelta[#[[1]],#\
[[2]]]&,Thread[{a,b}]]

GetIndexSets[\[Alpha]_List,a_List,b_List,y_List]:=
  Map[Flatten,
    {Position[(1-UnitStep[a-\[Alpha]])(1-UnitStep[-(b-\[Alpha])]),1],
      Position[MyKronDelta[\[Alpha],a]*UnitStep[y],1],
      Position[MyKronDelta[\[Alpha],b]*UnitStep[-y],1],
      Position[MyKronDelta[\[Alpha],b]*UnitStep[y],1],
      Position[MyKronDelta[\[Alpha],a]*UnitStep[-y],1]}]



UpperBoundarySet[Isets_]:=
  Isets[[1]] \[Union] Isets[[2]] \[Union] Isets[[3]]



LowerBoundarySet[Isets_]:=
  Isets[[1]] \[Union] Isets[[4]] \[Union] Isets[[5]]





ViolatingPairQ[Iup_,Ilow_,i_,j_,F_List,\[Tau]_]:=
  (MemberQ[Iup,i]\[And]
        MemberQ[Ilow,
          j]\[And](F[[j]]-
              F[[
                i]]>\[Tau]))\[Or](MemberQ[Ilow,i]\[And]
        MemberQ[Iup,
          j]\[And](F[[i]]-
              F[[j]]>\[Tau]))



GetViolatingPair[I0_,Iup_,Ilow_,F_List,\[Tau]_]:=
  Module[{o,i,j,l},
    If[Max[F[[I0]]]-
          Min[F[[I0]]]>\[Tau],
      o=Ordering[F[[I0]]];
      {I0[[First[o]]],
        I0[[Last[o]]]},
      (* else no violating pair on I0, scan entire {1...l} *)
      i=1;l=Length[F];
      While[i<=l,
        j=1;
        While[j<=l,
          If[ViolatingPairQ[Iup,Ilow,i,j,F,\[Tau]],
            Return[{i,j}]];
          j++];
        i++];
      Null]]





QPOptimalQ[Iup_,Ilow_,F_List,\[Tau]_]:=
  Max[F[[Ilow]]]-
      Min[F[[Iup]]]<=\[Tau]





RemainingOptions[r_List,opts___]:=
  Sequence[DeleteCases[{opts},Apply[Alternatives,Map[Rule[#,_]&,r]]]]



End[];

EndPackage[];
